#%%加载所需库文件
import numpy as np
import pandas as pd
from scipy import stats
np.set_printoptions(suppress=True)
#这个东西是让结果显示取消科学计数法，更加直观
def dottime(a,b):
      a=np.matrix(a)
      b=np.matrix(b)
      if np.size(a)!=np.size(b):
            print('传入矩阵的维度不同，无法进行数组运算！')
      else:
            c=np.zeros_like(a)
            c=np.matrix(c)
            for i in range(0,np.size(a,0)):
                  for j in range(0,np.size(a,1)):
                        c[i,j]=a[i,j]*b[i,j]
      return c
#和matlab里面的点乘操作一致！

#%%直接导入excle文件
excle_data=pd.read_excel(r"D:\my codes for matlab\self_learning_from_qingfeng\第14讲.主成分分析\代码和例题数据\作业数据.xlsx",sheet_name="Sheet1")
original_data=excle_data.drop(labels=['地区'],axis=1)
#该步骤代码的意思是删除原本表格中的第一列：地区，变为仅含有一个因变量和6个自变量的原始纯数表格

ori=original_data.drop(labels=['工业总产值y'],axis=1)

standard_data = stats.zscore(ori, axis=0,ddof=1)
#直接使用scipy中的stats.zscore就和matlab中的结果一样了
#之前我的代码应该就是数据的标准化出了问题
#问题在于在除以标准差的时候没有使用“减一”后的数据量作为计算标准差的个数
dmt_row=np.size(standard_data,axis=0)
dmt_col=np.size(standard_data,axis=1)

#第二步：计算标准化样本的协方差矩阵
covmat=np.cov(standard_data.T)
#np自带有协方差矩阵计算的函数
#但是要注意的是他默认每一行是一个特征，每一列是一个观测值，
#这与我们事先准备好的数据是相反的，所以要转置一下，.T

#第三步：计算特征值和特征向量
evalue, evictor = np.linalg.eig(covmat)
#特征值储存于evalue中，特征向量储存于evictor中

#print('特征值：{}\n特征向量：{}'.format(evalue, evictor))
#检查特征值、特征向量是否计算正确

#第四步：计算贡献率和累计贡献率

#贡献率实际上就是特征值占全体特征值的和的比例
dtb_pct=np.zeros_like(evalue)
for i in range(0,dmt_col):
      dtb_pct[i]=evalue[i]/np.sum(evalue)
#贡献率储存在dtb_pct中,distribution percentage
#下面计算累计贡献率：
adp=np.zeros_like(evalue)
asm=0
for i in range(0,dmt_col):
      asm += evalue[i]
      adp[i]= asm /np.sum(evalue)
#累计贡献率储存在adp中，accumulate distribution percentage
for i in range(0,dmt_col):
      print('第',i+1,'个指标对应的贡献率为：',dtb_pct[i])
      print('截止到该指标的累计贡献率为：',adp[i])

#%%第五步：提取主成分
#我们将主成分累计贡献率大于85%的前若干个特征值对应的指标定为提取出来的主成分。
#所以需要逐项查看累计贡献率中第几项开始超过0.85
k=int(input('请输入需要保留的主成分个数：'))
print('您选择保留了',k,'项主成分')
print('第',k,'项的贡献度为：',dtb_pct[k-1])
print('至此累计贡献度为：',adp[k-1])


F=np.zeros((dmt_col,k))   #储存所提取的主成分前面的系数
F=np.matrix(F)
for i in range(0,k):
      for j in range(0,dmt_col):
            F[j,i]  =  evictor[j,i]
print('您提取的主成分之前的系数矩阵已按照若干',k,'个列向量形式摆出如下：')
print(F)

#%%最后一步的计算，计算主成分的值
G=np.ones([dmt_row,k])
G=np.matrix(G)#用于计算最终的主成分结果



for l in range(0,k):
      ai=F[:,l].T
      Ai=np.tile(ai,(dmt_row,1))
      G[:,l]=np.sum((dottime(Ai,standard_data)),1)

#注！np提供了np.tile的函数，与matlab里面的repmat作用完全相同！

results=pd.DataFrame(G)
#写入数据集，便于后续导出至excle中

# %%保存导出数据
writer = pd.ExcelWriter('~\desktop\作业求得的主成分.xlsx')
#新建一个writer，用于设定导出的excle文件的位置
results.to_excel(writer,'主成分取值')
F=pd.DataFrame(F)
F.to_excel(writer,'主成分系数')
#导出文件
writer.save()
#保存文件
print('您所需要的东西已经导出至桌面，请查收')


# %%
